library(lme4)
data = read.csv('human_computer_survey_data.csv',header=TRUE)

# --------------------
# Remove some data
# --------------------
late_learners = data[data$learn_english>3,]                          # Late learners of English (age > 3; these are one hindi and one chinese speaker)
short_pres    = data[(data$read_time/data$doc_length)<10,]           # Responses on texts that were visible for less than 10ms per character
topic_error   = data[data$subj_answer_correct==0,]                   # Tesponses on texts the where topic question was not answered correctly

data = data[setdiff(rownames(data),rownames(late_learners)),]
data = data[setdiff(rownames(data),rownames(short_pres)),]
data = data[setdiff(rownames(data),rownames(topic_error)),]

# --------------------------------------------
# Select data subsets and recode binary values
# --------------------------------------------

data$doc_ordinal        = data$doc_ordinal-1                    # Count doc_ordinal from 0, so main effect of response_type/response_system will be the estimate before any learning has taken place
data$response_type      = 2*data$response_type-1                # Human: -1; Coputer: +1
humans                  = data[data$response_type==-1,]          # Human-generated responses only
machine                 = data[data$response_type==1,]          # Computer-generated responses only
machine$response_system = 2*machine$response_system-3           # Without KB: -1; With KB: +1  
KB                      = machine[machine$response_system==1,]
noKB                    = machine[machine$response_system==-1,]

#-----------------------------------
# Compute score averages and 95% CIs
#-----------------------------------

human_mean_CI = c(mean(humans$response_score),  1.96*sd(humans$response_score) /sqrt(nrow(humans)))
comp_mean_CI  = c(mean(machine$response_score), 1.96*sd(machine$response_score)/sqrt(nrow(machine)))
KB_mean_CI    = c(mean(KB$response_score),      1.96*sd(KB$response_score)     /sqrt(nrow(KB)))
noKB_mean_CI  = c(mean(noKB$response_score),    1.96*sd(noKB$response_score)   /sqrt(nrow(noKB)))

#--------------------------------------------------------------------------
# Fit regression models: human versus machine, and system 1 versus system 2
#--------------------------------------------------------------------------

# Include by-participant random slopes of response_type (some participants may be more sensitive to the difference than others)
# and of doc_ordinal (some participants may learn faster than others).
model_human_machine = summary(lmer(response_score ~ response_type*doc_ordinal + (1|assignment_id) + (0+response_type|assignment_id) + (0+doc_ordinal|assignment_id), data))

# Include by-participant random slopes of response_system, response_sentiment, abs(response_sentiment) and doc_ordinal
model_machines      = summary(lmer(response_score ~ response_system*doc_ordinal + response_sentiment + abs(response_sentiment) + (1|assignment_id) + (0+response_system|assignment_id) + (0+doc_ordinal|assignment_id) + (0+response_sentiment|assignment_id) + (0+abs(response_sentiment)|assignment_id), machine))

#-----------------------------------------------------------------
# Do MCMC sampling to get posterior distribution over coefficients
#-----------------------------------------------------------------

Pb0_human_machine = array(dim=c(1,3))
set.seed(1)
mcmc_human_machine = mcmcsamp(model_human_machine, n=10000)
HPD_human_machine  = HPDinterval(mcmc_human_machine)
Pb0_human_machine[1] = mean(mcmc_human_machine@fixef[2,]<=0)
Pb0_human_machine[2] = mean(mcmc_human_machine@fixef[3,]<=0)
Pb0_human_machine[3] = mean(mcmc_human_machine@fixef[4,]<=0)

Pb0_machines = array(dim=c(1,5))
set.seed(1)
mcmc_machines = mcmcsamp(model_machines, n=10000)
HPD_machines  = HPDinterval(mcmc_machines)
Pb0_machines[1] = mean(mcmc_machines@fixef[2,]<=0)
Pb0_machines[2] = mean(mcmc_machines@fixef[3,]<=0)
Pb0_machines[3] = mean(mcmc_machines@fixef[4,]<=0)
Pb0_machines[4] = mean(mcmc_machines@fixef[5,]<=0)
Pb0_machines[5] = mean(mcmc_machines@fixef[6,]<=0)

save.image('analysis_Exp1.RData')